clear;

listDir = dir('Data*');
saveArg = 1;
numBins = 72;
lowerThreshold = 10^-5;
intrudeErupt = zeros(15,1);
time = 	([7000:3500:52500,55600] - 3262)/10;


%% Processing

L = nan(length(listDir),1);
B = L;
thetaMean = L;
thetaMedian = L;
thetaDev = L;
thetaMad = L;
correctedMean = L;
correctedMedian = L;
correctedDev = L;
correctedMad = L;
reducedMean = L;
reducedMedian = L;
reducedDev = L;
reducedMad = L;
reducedCorMean = L;
reducedCorMedian = L;
reducedCorDev = L;
reducedCorMad = L;
vectorDistMed = nan(numBins,length(listDir));
vectorDistMad = nan(numBins,length(listDir));
correctedDistMed = nan(numBins,length(listDir));
correctedDistMad = nan(numBins,length(listDir));
reducedDistMed = nan(numBins,length(listDir));
reducedDistMad = nan(numBins,length(listDir));
velocityMean = L;

for i = 1:length(listDir)
    data = importdata(sprintf('%s\\B0001.dat',listDir(i).name));
    data = data.data;
    x = data(:,1);
    y = data(:,2);
    vx = data(:,3);
    vy = data(:,4);
    noInd = find(vx==0&vy==0);
    x(noInd) = [];
    y(noInd) = [];
    vx(noInd) = [];
    vy(noInd) = [];
    L(i) = max(y) - min(y);
    B(i) = max(x) - min(x);
    
    vector = sqrt(vx.^2+vy.^2);
    theta = acos(vx./vector);
    
    correctedTheta = theta;
    correctedTheta(vy<0) = 2*pi - correctedTheta(vy<0);
    correctedTheta = correctedTheta + pi/2;
    correctedTheta(correctedTheta>2*pi) = correctedTheta(correctedTheta>2*pi) - 2*pi;
    %     correctedTheta = pi - theta(vx<0);
    
%     lowerThreshold = max(vector)/100*5;
    threshold = max([lowerThreshold,mean(vector)-std(vector)]);
%     slowInd = find(vector<threshold|vector>1e-4);
%     reducedTheta = theta;
%     reducedCorrected = correctedTheta;
%     reducedVector = vector;
%     reducedVy = vy;
%     reducedTheta(slowInd) = [];
%     reducedCorrected(slowInd) = [];
%     reducedVector(slowInd) = [];
%     reducedVy(slowInd) = [];
    
    thetaMean(i) = mean(theta);
    thetaMedian(i) = median(theta);
    thetaDev(i) = std(theta);
    thetaMad(i) = mad(theta);
    correctedMean(i) = mean(correctedTheta);
    correctedMedian(i) = median(correctedTheta);
    correctedDev(i) = std(correctedTheta);
    correctedMad(i) = mad(correctedTheta);
%     reducedMean(i) = mean(reducedTheta);
%     reducedMedian(i) = median(reducedTheta);
%     reducedDev(i) = std(reducedTheta);
%     reducedMad(i) = mad(reducedTheta);
%     reducedCorMean(i) = mean(reducedCorrected);
%     reducedCorMedian(i) = median(reducedCorrected);
%     reducedCorDev(i) = std(reducedCorrected);
%     reducedCorMad(i) = mad(reducedCorrected);
    
    flagInd = ones(size(theta));
    for j = 2:numBins
        checkInd = find(theta>=2*pi*(j-1)/numBins & ...
           theta<=2*pi*j/numBins);
        vectorDistMed(j,i) = median(vector(checkInd));
        vectorDistMad(j,i) = mad(vector(checkInd));

        checkInd = find(correctedTheta>=2*pi*(j-1)/numBins & ...
           correctedTheta<=2*pi*j/numBins);
        correctedDistMed(j,i) = median(vector(checkInd));
        correctedDistMad(j,i) = mad(vector(checkInd));

%         checkInd = find(reducedTheta>=2*pi*(j-1)/numBins & ...
%            reducedTheta<=2*pi*j/numBins);
%         reducedDistMed(j,i) = median(reducedVector(checkInd));
%         reducedDistMad(j,i) = mad(reducedVector(checkInd));
        
        checkInd = find(theta>=2*pi*(j-1)/numBins & ...
           theta<=2*pi*j/numBins & ...
           vector<vectorDistMed(j,i)*4 & vector>vectorDistMed(j,i)/4 & ...
           vector>threshold);
        if isempty(checkInd) == 0
            flagInd(checkInd) = 0;
        end
    end
    thetaBins = [2*pi/numBins:2*pi/numBins:2*pi] - 2*pi/numBins/2;
    
    reducedTheta = theta(flagInd==0);
    reducedVector = vector(flagInd==0);
    reducedCorrected = correctedTheta(flagInd==0);
    reducedVy = vy(flagInd==0);
    for j = 2:numBins
        checkInd = find(reducedTheta>=2*pi*(j-1)/numBins & ...
           reducedTheta<=2*pi*j/numBins);
        reducedDistMed(j,i) = median(reducedVector(checkInd));
        reducedDistMad(j,i) = mad(reducedVector(checkInd));
    end
    velocityMean(i) = mean(reducedVector);
    reducedMean(i) = mean(reducedTheta);
    reducedMedian(i) = median(reducedTheta);
    reducedDev(i) = std(reducedTheta);
    reducedMad(i) = mad(reducedTheta);
    reducedCorMean(i) = mean(reducedCorrected);
    reducedCorMedian(i) = median(reducedCorrected);
    reducedCorDev(i) = std(reducedCorrected);
    reducedCorMad(i) = mad(reducedCorrected);
    
    
    figure(1); clf;
    set(gcf,'Units','normalized','OuterPosition',[0.15 0.05 0.45 0.95]);
    subplot(2,2,1);
        rose(theta,numBins);
    subplot(2,2,2);
        rose(correctedTheta-pi/2,numBins);
    subplot(2,2,3);
        rose(reducedTheta,numBins);
    subplot(2,2,4);
        rose(reducedCorrected-pi/2,numBins);
    if saveArg == 1
        print(gcf,sprintf('Rose_diagrams\\Rose_%s',listDir(i).name),...
            '-djpeg','-r300');
    end
    
    figure(2); clf;
    set(gcf,'Units','normalized','OuterPosition',[0.1 0.07 0.7 0.9]);
    subplot(2,2,1);
        plot(theta*180/pi,vector,'b.'); hold on;
        errorbar(thetaBins*180/pi,vectorDistMed(:,i),vectorDistMad(:,i),...
            'r.-','Linewidth',2);
    subplot(2,2,2);
        plot(correctedTheta*180/pi,vector,'b.'); hold on;
        errorbar(thetaBins*180/pi,correctedDistMed(:,i),correctedDistMad(:,i),...
            'r.-','Linewidth',2);
    subplot(2,2,3);
        plot(reducedTheta*180/pi,reducedVector,'b.'); hold on;
        errorbar(thetaBins*180/pi,reducedDistMed(:,i),reducedDistMad(:,i),...
            'r.-','Linewidth',2);
    subplot(2,2,4);
        plot(thetaBins*180/pi,vectorDistMed(:,i),...
            'b.-','Linewidth',2); hold on;
        plot(thetaBins*180/pi,correctedDistMed(:,i),...
            'g.-','Linewidth',2);
        plot(thetaBins*180/pi,reducedDistMed(:,i),...
            'r.-','Linewidth',2);
        set(gca,'YScale','log');
    for j = 1:4
        subplot(2,2,j);
        ylabel('v (m/s)'); xlabel('\theta (deg)');
        if j <= 3
            ylim([0,1e-4]);
        end
    end
    if saveArg == 1
        print(gcf,sprintf('Rose_diagrams\\ThetaChart_%s',listDir(i).name),...
            '-djpeg','-r300');
    end
    
    figure(3); clf;
    set(gcf,'Units','normalized','OuterPosition',[0.15 0.1 0.65 0.85]);
    subplot(2,2,1);
        hist(theta,72); hold on;
        plot(thetaMean(i)*[1,1],get(gca,'YLim'),'k-');
        plot((thetaMean(i)-thetaDev(i))*[1,1],get(gca,'YLim'),'k--');
        plot((thetaMean(i)+thetaDev(i))*[1,1],get(gca,'YLim'),'k--');
        plot(thetaMedian(i)*[1,1],get(gca,'YLim'),'r-');
        plot((thetaMedian(i)-thetaMad(i))*[1,1],get(gca,'YLim'),'r--');
        plot((thetaMedian(i)+thetaMad(i))*[1,1],get(gca,'YLim'),'r--');
    subplot(2,2,2);
        hist(correctedTheta,72); hold on;
        plot(correctedMean(i)*[1,1],get(gca,'YLim'),'k-');
        plot((correctedMean(i)-correctedDev(i))*[1,1],get(gca,'YLim'),'k--');
        plot((correctedMean(i)+correctedDev(i))*[1,1],get(gca,'YLim'),'k--');
        plot(correctedMedian(i)*[1,1],get(gca,'YLim'),'r-');
        plot((correctedMedian(i)-correctedMad(i))*[1,1],get(gca,'YLim'),'r--');
        plot((correctedMedian(i)+correctedMad(i))*[1,1],get(gca,'YLim'),'r--');
    subplot(2,2,3);
        hist(reducedTheta,72); hold on;
        plot(reducedMean(i)*[1,1],get(gca,'YLim'),'k-');
        plot((reducedMean(i)-reducedDev(i))*[1,1],get(gca,'YLim'),'k--');
        plot((reducedMean(i)+reducedDev(i))*[1,1],get(gca,'YLim'),'k--');
        plot(reducedMedian(i)*[1,1],get(gca,'YLim'),'r-');
        plot((reducedMedian(i)-reducedMad(i))*[1,1],get(gca,'YLim'),'r--');
        plot((reducedMedian(i)+reducedMad(i))*[1,1],get(gca,'YLim'),'r--');
    subplot(2,2,4);
        hist(reducedCorrected,72); hold on;
        plot(reducedCorMean(i)*[1,1],get(gca,'YLim'),'k-');
        plot((reducedCorMean(i)-reducedCorDev(i))*[1,1],get(gca,'YLim'),'k--');
        plot((reducedCorMean(i)+reducedCorDev(i))*[1,1],get(gca,'YLim'),'k--');
        plot(reducedCorMedian(i)*[1,1],get(gca,'YLim'),'r-');
        plot((reducedCorMedian(i)-reducedCorMad(i))*[1,1],get(gca,'YLim'),'r--');
        plot((reducedCorMedian(i)+reducedCorMad(i))*[1,1],get(gca,'YLim'),'r--');
    if saveArg == 1
        print(gcf,sprintf('Rose_diagrams\\ThetaHist_%s',listDir(i).name),...
            '-djpeg','-r300');
    end
    
    if i == 1 || i == 11
        figure(4);
        if i == 1
            set(gcf,'units','normalized','OuterPosition',[0.16 0.13 0.50 0.77]);
            subplot(2,2,1);
        elseif i == 11
            subplot(2,2,3);
        end
        scatter(x,y,10,vx*1000,'filled'); c=colorbar;
        set(gca,'XTick',[],'YTick',[]);
        ylabel(c,'v_x (mm/s)');
        caxis(max(vx)*1000*[-1,1]);
        colormap jet
        if i == 1
            subplot(2,2,2);
        elseif i == 11
            subplot(2,2,4);
        end
        hist(reducedTheta,72); hold on;
        plot(reducedMean(i)*[1,1],get(gca,'YLim'),'k-');
        plot((reducedMean(i)-reducedDev(i))*[1,1],get(gca,'YLim'),'k--');
        plot((reducedMean(i)+reducedDev(i))*[1,1],get(gca,'YLim'),'k--');
        plot(reducedMedian(i)*[1,1],get(gca,'YLim'),'r-');
        plot((reducedMedian(i)-reducedMad(i))*[1,1],get(gca,'YLim'),'r--');
        plot((reducedMedian(i)+reducedMad(i))*[1,1],get(gca,'YLim'),'r--');
        xlim([0,pi]);
        set(gca,'Xtick',[0,pi/2,pi],'XTickLabel',[0,90,180]);
            
        if i == 11 && saveArg == 1
            print(gcf,'Flow_samples','-djpeg','-r300');
            saveas(gcf,'Flow_samples','fig');
        end
    end
        
    
    pause(0.1);
end


%% Output

figure(10); clf;
set(gcf,'Units','normalized','OuterPosition',[0.05 0.15 0.85 0.55]);
subplot(1,3,1);
    errorbar(thetaMean*180/pi,thetaDev*180/pi,'b.-','Markersize',10); hold on;
    errorbar(correctedMean*180/pi,correctedDev*180/pi,'g.-','Markersize',10);
    errorbar(reducedMean*180/pi,reducedDev*180/pi,'r.-','Markersize',10);
    errorbar(reducedCorMean*180/pi,reducedCorDev*180/pi,'k.-','Markersize',10);
    errorbar(thetaMedian*180/pi,thetaMad*180/pi,'b.--','Markersize',10); hold on;
    errorbar(correctedMedian*180/pi,correctedMad*180/pi,'g.--','Markersize',10);
    errorbar(reducedMedian*180/pi,reducedMad*180/pi,'r.--','Markersize',10);
    errorbar(reducedCorMedian*180/pi,reducedCorMad*180/pi,'k.--','Markersize',10);
    xlabel('Frame #'); ylabel('\theta (deg)'); title('Mean, Deviation');
subplot(1,3,2);
    plot(thetaMean*180/pi,'b.-','Markersize',10); hold on;
    plot(correctedMean*180/pi,'g.-','Markersize',10);
    plot(reducedMean*180/pi,'r.-','Markersize',10);
    plot(reducedCorMean*180/pi,'k.-','Markersize',10);
    plot(thetaMedian*180/pi,'b.--','Markersize',10);
    plot(correctedMedian*180/pi,'g.--','Markersize',10);
    plot(reducedMedian*180/pi,'r.--','Markersize',10);
    plot(reducedCorMedian*180/pi,'k.--','Markersize',10);
    xlabel('Frame #'); title('Mean');
subplot(1,3,3);
    plot(thetaDev*180/pi,'b.-','Markersize',10); hold on;
    plot(correctedDev*180/pi,'g.-','Markersize',10);
    plot(reducedDev*180/pi,'r.-','Markersize',10);
    plot(reducedCorDev*180/pi,'k.-','Markersize',10);
    plot(thetaMad*180/pi,'b.--','Markersize',10);
    plot(correctedMad*180/pi,'g.--','Markersize',10);
    plot(reducedMad*180/pi,'r.--','Markersize',10);
    plot(reducedCorMad*180/pi,'k.--','Markersize',10);
    xlabel('Frame #'); title('Deviation');
    if saveArg == 1
        print(gcf,'Theta_Dist','-djpeg','-r300');
    end

if saveArg == 1
    save('Output.mat','time','B','L','x','y','vx','vy',...
        'corrected*','reduced*','theta*','vector*','intrudeErupt','velocityMean');
end



